
% ZM: 泽尼克多项式矩阵；mask：采样点阵；D：泽尼克圆的直径；S：mask的面积
function M = OrthogonalMatrix(ZM, D, mask, S)

lnstp = D/size(ZM,1); %采样间隔
%NN = size(ZM,1)*size(ZM,2);
j_max = size(ZM,3);
%NS = length(find(mask));
%S = NS/(NN-2*sqrt(NN)) * D^2;

CZZ = zeros(j_max,j_max);
for j = 1:j_max
    Z1 = ZM(:,:,j);
    for k = 1:j
        Z2 = ZM(:,:,k);
       if j == k
           CZZ(j,k) = InnerProduct(Z1, Z2);
       else
           CZZ(j,k) = InnerProduct(Z1, Z2);
           CZZ(k,j) = CZZ(j,k);
       end
    end
end

Q = chol(CZZ,'upper'); % Q'Q=CZZ
% [U,V] = eig(CZZ);
% Q = real(U*sqrt(V));
M = inv(Q');

    % 定义一个子函数，两个ZM函数在mask区域内做内积
    function ZZ = InnerProduct(Z1, Z2)
        ZZ = sum(sum(Z1(mask).*Z2(mask)))*lnstp^2/S;
    end

% end of the function
end